Graphs

Load Estimates

wv4_est_weighted <- readRDS(here("models", "Wave 4", 
                                 "Weighted Bayesian", "wv4_weighted_est.RDS"))
wv7_est_weighted <- readRDS(here("models", "Wave 7", 
                                 "Weighted Bayesian", "wv7_weighted_est.RDS"))
obs_wv4_weighted <- readRDS(here("models", "Wave 4", 
                                 "Weighted Bayesian", "obs_wv4_weighted.RDS"))
obs_wv7_weighted <- readRDS(here("models", "Wave 7", 
                                 "Weighted Bayesian", "obs_wv7_weighted.RDS"))

Data Prep

# Wave 4
wv4_est_weighted <- left_join(wv4_est_weighted, obs_wv4_weighted, by = "strata")

wv4_est_weighted <- wv4_est_weighted %>% 
 mutate(wave = 4) %>% 
 mutate(est_obs_dif = eprod_edsd - mean) %>% 
 mutate(mean = mean*100,
        q5 = q5*100,
        q95 = q95*100,
        eprod_edsd = eprod_edsd*100,
        est_obs_dif = est_obs_dif*100)

# Wave 7
wv7_est_weighted <- left_join(wv7_est_weighted, obs_wv7_weighted, by = "strata")

wv7_est_weighted <- wv7_est_weighted %>% 
 mutate(wave = 7) %>% 
 mutate(est_obs_dif = eprod_edsd - mean) %>% 
 mutate(mean = mean*100,
        q5 = q5*100,
        q95 = q95*100,
        eprod_edsd = eprod_edsd*100,
        est_obs_dif = est_obs_dif*100)

# Dataframe with both waves, rows
est_weighted <- rbind.data.frame(wv4_est_weighted, wv7_est_weighted)

# Dataframe with both waves, columns
est_weighted_wide <- est_weighted %>%
  pivot_wider(
    id_cols = c(strata, strata_label),
    names_from = wave,
    values_from = c(mean, q5, q95, eprod_edsd),
    names_glue = "{.value}_wave{wave}") %>% 
 mutate(est_wave_diff = mean_wave7-mean_wave4)

Wave 4: Weighted Estimates

Etimates graph

wv4_est_weighted %>%
 ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum Rank") +
  coord_flip() +
  theme_bw()

Faceted by region

wv4_est_weighted %>%
  mutate(region = word(strata_label, 1, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ region, scales = "free_y", ncol = 1)

Faceted by race/ethnicity

wv4_est_weighted %>%
  mutate(race_ethnicity = word(strata_label, 2, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ race_ethnicity, scales = "free_y", ncol = 1)

Faceted by sex

wv4_est_weighted %>%
  mutate(sex = word(strata_label, 3, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ sex, scales = "free_y", ncol = 1)

Faceted by age

wv4_est_weighted %>%
  mutate(age = word(strata_label, 4, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ age, scales = "free_y", ncol = 1)

Estimated vs. observed graph

On average, the mean predicted EDSD rate was only 0.14 percentage points low than the weighted, observed rate (5.44% vs. 5.47%).

wv4_est_weighted %>%
  ggplot(aes(x = reorder(strata_label, abs(est_obs_dif)))) +
  geom_point(aes(y = mean, color = "Predicted"), size = 2) +
  geom_point(aes(y = eprod_edsd, color = "Observed"), size = 2) +
  scale_color_manual(values = c("Predicted" = "black", "Observed" = "purple")) +
  ylab("Percent EDSD E-Product Use") +
  xlab("Stratum (Ranked by Accuracy, Lowest to Highest)") +
  coord_flip() +
  theme_bw() +
  theme(legend.title = element_blank())

describe(wv4_est_weighted$eprod_edsd)

##    vars  n mean   sd median trimmed  mad min   max range skew kurtosis   se
## X1    1 96 5.47 4.37   4.45    4.79 3.52   0 22.02 22.02 1.44     1.94 0.45

describe(wv4_est_weighted$mean)

##    vars  n mean   sd median trimmed  mad  min  max range skew kurtosis   se
## X1    1 96 5.33 3.58   4.89    4.91 3.94 1.25 16.4 15.15 1.06     0.85 0.37

Wave 7: Weighted Estimates

Full

wv7_est_weighted %>%
 ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted Percent EDSD E-Product User") +
  xlab("Stratum Rank") +
  coord_flip() +
  theme_bw()

Faceted by region

wv7_est_weighted %>%
  mutate(region = word(strata_label, 1, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted Percent EDSD E-Product User") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ region, scales = "free_y", ncol = 1)

Faceted by race/ethnicity

wv7_est_weighted %>%
  mutate(race_ethnicity = word(strata_label, 2, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted Percent EDSD E-Product User") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ race_ethnicity, scales = "free_y", ncol = 1)

Faceted by sex

wv7_est_weighted %>%
  mutate(sex = word(strata_label, 3, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted Percent EDSD E-Product User") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ sex, scales = "free_y", ncol = 1)

Faceted by age

wv7_est_weighted %>%
  mutate(age = word(strata_label, 4, sep = fixed(","))) %>% 
  ggplot(aes(y = mean, x = reorder(strata_label, mean))) +
  geom_pointrange(aes(ymin = q5, ymax = q95)) +
  ylab("Predicted Percent EDSD E-Product User") +
  xlab("Stratum") +
  coord_flip() +
  theme_bw() +
  facet_wrap(~ age, scales = "free_y", ncol = 1)

Estimated vs. observed graph

On average, the mean predicted EDSD rate was only 0.17 percentage points higher than the weighted, observed rate (9.67% vs. 9.5%).

wv7_est_weighted %>%
  ggplot(aes(x = reorder(strata_label, abs(est_obs_dif)))) +
  geom_point(aes(y = mean, color = "Predicted"), size = 2) +
  geom_point(aes(y = eprod_edsd, color = "Observed"), size = 2) +
  scale_color_manual(values = c("Predicted" = "black", "Observed" = "purple")) +
  ylab("Predicted EDSD E-Product Use Rates") +
  xlab("Stratum (Ranked by Accuracy, Lowest to Highest)") +
  coord_flip() +
  theme_bw() +
  theme(legend.title = element_blank())

describe(wv7_est_weighted$eprod_edsd)

##    vars  n mean  sd median trimmed  mad min   max range skew kurtosis   se
## X1    1 96 9.67 7.1   8.61    9.14 8.19   0 24.71 24.71 0.53     -0.9 0.72

describe(wv7_est_weighted$mean)

##    vars  n mean   sd median trimmed mad  min   max range skew kurtosis   se
## X1    1 96  9.5 6.61   8.98    9.01 9.6 1.19 23.65 22.46 0.37       -1 0.67

Wave 4 vs. Wave 7 Estimates

est_weighted_wide %>%
  ggplot(aes(x = reorder(strata_label, est_wave_diff))) +
  
  # Wave 4 point + interval
  geom_point(aes(y = mean_wave4, color = "Wave 4"), size = 2) +
  geom_pointrange(aes(y = mean_wave4, ymin = q5_wave4, ymax = q95_wave4), color = "turquoise4") +

  # Wave 7 point + interval
  geom_point(aes(y = mean_wave7, color = "Wave 7"), size = 2) +
  geom_pointrange(aes(y = mean_wave7, ymin = q5_wave7, ymax = q95_wave7), color = "darkorchid4") +

  scale_color_manual(values = c("Wave 4" = "turquoise4", "Wave 7" = "darkorchid4")) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum (Ranked by Change, Greatest to Least)") +
  coord_flip() +
  theme_bw() +
  theme(legend.title = element_blank())

Faceted by region

region_facet <- est_weighted_wide %>%
  mutate(region = word(strata_label, 1, sep = fixed(","))) %>% 
  ggplot(aes(x = reorder(strata_label, est_wave_diff))) +
  geom_point(aes(y = mean_wave4, color = "Wave 4"), size = 2) +
  geom_pointrange(aes(y = mean_wave4, ymin = q5_wave4, ymax = q95_wave4), color = "turquoise4") +
  geom_point(aes(y = mean_wave7, color = "Wave 7"), size = 2) +
  geom_pointrange(aes(y = mean_wave7, ymin = q5_wave7, ymax = q95_wave7), color = "darkorchid4") +
  scale_color_manual(values = c("Wave 4" = "turquoise4", "Wave 7" = "darkorchid4")) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum (Ranked by Change, Greatest to Least)") +
  coord_flip() +
  facet_wrap(~ region, scales = "free_y", ncol = 1) +
  theme_bw() +
  theme(legend.title = element_blank())

region_facet

Faceted by race/ethnicity

est_weighted_wide %>%
  mutate(race_ethnicity = word(strata_label, 2, sep = fixed(","))) %>% 
  ggplot(aes(x = reorder(strata_label, est_wave_diff))) +
  geom_point(aes(y = mean_wave4, color = "Wave 4"), size = 2) +
  geom_pointrange(aes(y = mean_wave4, ymin = q5_wave4, ymax = q95_wave4), color = "turquoise4") +
  geom_point(aes(y = mean_wave7, color = "Wave 7"), size = 2) +
  geom_pointrange(aes(y = mean_wave7, ymin = q5_wave7, ymax = q95_wave7), color = "darkorchid4") +
  scale_color_manual(values = c("Wave 4" = "turquoise4", "Wave 7" = "darkorchid4")) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum (Ranked by Change, Greatest to Least)") +
  coord_flip() +
  facet_wrap(~ race_ethnicity, scales = "free_y", ncol = 1) +
  theme_bw() +
  theme(legend.title = element_blank())

Faceted by sex

est_weighted_wide %>%
  mutate(sex = word(strata_label, 3, sep = fixed(","))) %>% 
  ggplot(aes(x = reorder(strata_label, est_wave_diff))) +
  geom_point(aes(y = mean_wave4, color = "Wave 4"), size = 2) +
  geom_pointrange(aes(y = mean_wave4, ymin = q5_wave4, ymax = q95_wave4), color = "turquoise4") +
  geom_point(aes(y = mean_wave7, color = "Wave 7"), size = 2) +
  geom_pointrange(aes(y = mean_wave7, ymin = q5_wave7, ymax = q95_wave7), color = "darkorchid4") +
  scale_color_manual(values = c("Wave 4" = "turquoise4", "Wave 7" = "darkorchid4")) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum (Ranked by Change, Greatest to Least)") +
  coord_flip() +
  facet_wrap(~ sex, scales = "free_y", ncol = 1) +
  theme_bw() +
  theme(legend.title = element_blank())

Faceted by age

age_facet <- est_weighted_wide %>%
  mutate(age = word(strata_label, 4, sep = fixed(","))) %>% 
  ggplot(aes(x = reorder(strata_label, est_wave_diff))) +
  geom_point(aes(y = mean_wave4, color = "Wave 4"), size = 2) +
  geom_pointrange(aes(y = mean_wave4, ymin = q5_wave4, ymax = q95_wave4), color = "turquoise4") +
  geom_point(aes(y = mean_wave7, color = "Wave 7"), size = 2) +
  geom_pointrange(aes(y = mean_wave7, ymin = q5_wave7, ymax = q95_wave7), color = "darkorchid4") +
  scale_color_manual(values = c("Wave 4" = "turquoise4", "Wave 7" = "darkorchid4")) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum (Ranked by Change, Greatest to Least)") +
  coord_flip() +
  facet_wrap(~ age, scales = "free_y", ncol = 1) +
  theme_bw() +
  theme(legend.title = element_blank())

age_facet

Faceted by Age and Sex

age_sex_facet <- est_weighted_wide %>%
  mutate(age = word(strata_label, 4, sep = fixed(",")),
         sex = word(strata_label, 3, sep = fixed(","))) %>% 
  ggplot(aes(x = reorder(strata_label, est_wave_diff))) +
  geom_point(aes(y = mean_wave4, color = "Wave 4"), size = 2) +
  geom_pointrange(aes(y = mean_wave4, ymin = q5_wave4, ymax = q95_wave4), color = "turquoise4") +
  geom_point(aes(y = mean_wave7, color = "Wave 7"), size = 2) +
  geom_pointrange(aes(y = mean_wave7, ymin = q5_wave7, ymax = q95_wave7), color = "darkorchid4") +
  scale_color_manual(values = c("Wave 4" = "turquoise4", "Wave 7" = "darkorchid4")) +
  ylab("Predicted EDSD E-Product Use (95% Crl)") +
  xlab("Stratum (Ranked by Change, Greatest to Least)") +
  coord_flip() +
  facet_wrap(~ interaction(age, sex), scales = "free_y") +
  theme_bw() +
  theme(legend.title = element_blank())

age_sex_facet